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ABSTRACT 


The engulfment of substellar bodies (SBs, such as brown dwarfs and planets) by giant stars is a possible 
explanation for rapidly rotating giants, lithium-rich giants, and the presence of SBs in close orbits around 
subdwarfs and white dwarfs. We simulate the flow in the vicinity of an engulfed SB in three-dimensional 
hydrodynamics. We model the SB as a rigid body with a reflective surface because it cannot accrete. This 
reflective boundary changes the flow morphology to resemble that of engulfed compact objects with outflows. 
We measure the drag coefficients for the ram pressure and gravitational drag forces acting on the SB, and use 
them to integrate its trajectory inside the star. We find that engulfment can increase the luminosity of a 1M 
star by up to a few orders of magnitude. The time for the star to return to its original luminosity is up to a few 
thousand years when the star has evolved to ~ 10; and up to a few decades at the tip of the red giant branch. 
No SBs can eject the envelope of a 1 Mo star before it evolves to + 10 R5, if the orbit of the SB is the only 
energy source contributing to the ejection. In contrast, SBs as small as ~ 10Mjup can eject the envelope at the tip 
of the red giant branch. The numerical framework we introduce here can be used to study planetary engulfment 


in a simplified setting that captures the physics of the flow at the scale of the SB. 


1. INTRODUCTION 


Common-envelope evolution (hereafter CEE; Paczynski 
1976) is a process in which a star engulfs a companion (sub- 
stellar or otherwise). The known planetary system architec- 
tures imply that a large fraction of planets and brown dwarfs 
(hereafter substellar bodies, SBs) will eventually undergo CEE 
(Villaver & Livio 2009; Mustill & Villaver 2012; Nordhaus 
& Spiegel 2013; Schlaufman & Winn 2013; Sun et al. 2018). 
Throughout this work, we will refer to CEE between a star 
and an SB as “planetary engulfment,” and use “CEE” for the 
more general interaction between a star and a companion of 
any mass. 

Planetary engulfment is a possible explanation for several 
unsolved problems in stellar and planetary system evolu- 
tion. Observations have found SBs in close orbits around 
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subdwarfs and white dwarfs (Schmidt et al. 2005; Littlefair 
et al. 2006; Maxted et al. 2006; Littlefair et al. 2007; Sil- 
vestri et al. 2007; Littlefair et al. 2008; Geier et al. 2009; 
Charpinet et al. 2011; Breedt et al. 2012; Casewell et al. 
2012; Liu et al. 2012; Rebassa-Mansergas et al. 2012; Beuer- 
mann et al. 2013; Steele et al. 2013; McAllister et al. 2015; 
Almeida et al. 2017; Schaffenroth et al. 2015; Parsons et al. 
2017; Pala et al. 2018; Casewell et al. 2020; Vanderburg et al. 
2020; Schaffenroth et al. 2021; van Roestel et al. 2021, for 
a summary see Kruckow et al. 2021). These systems might 
have reached their current orbital configurations dynamically 
through the Kozai-Lidov mechanism (Kozai 1962; Fabrycky 
& Tremaine 2007; Katz et al. 2011; Naoz et al. 2012; Socrates 
et al. 2012; Shappee & Thompson 2013; Mufioz & Petrovich 
2020; O'Connor et al. 2021) or via an engulfment phase in 
which the SB ejected the envelope of the star that engulfed it 
(Livio & Soker 1984; Nelemans & Tauris 1998; Lagos et al. 
2021; Merlov et al. 2021; Zorotovic & Schreiber 2022). Dur- 
ing engulfment, orbital energy dissipation shrinks the orbit of 
the system significantly. Even if the SB does not survive, en- 
gulfment might result in an isolated white dwarf with > MG 
magnetic fields (Nordhaus et al. 2011; Guidarelli et al. 2019). 
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Engulfment could explain observations of anomalous ro- 
tation among some giant stars. During engulfment, the SB 
transfers the angular momentum of its orbit into the stellar 
envelope, resulting in either enhanced or reduced rotation, 
depending on the alignment of the angular momentum vec- 
tors of the star and the orbit. Engulfment can speed up the 
surface of giant stars up to the observed values, and even up 
to a significant fraction of their critical speeds (Peterson et al. 
1983; Soker 1998; Siess & Livio 1999; Zhang & Penev 2014; 
Privitera et al. 2016a,b; Qureshi et al. 2018; Stephan et al. 
2020). 

While the stellar surface abundance of the "Li isotope gen- 
erally decreases throughout stellar evolution (Bodenheimer 
1965; Deliyannis et al. 2000; Piau & Turck-Chiéze 2002; Bau- 
mann et al. 2010; Monroe et al. 2013; Meléndez et al. 2014; 
Carlos et al. 2016, 2019; Soares-Furtado et al. 2021), drop- 
ping significantly at the onset of the first dredge-up phase, 
~ 1% of giants have abundances > 1.5 dex (e.g., Wallerstein 
& Sneden 1982; Brown et al. 1989; Balachandran et al. 2000; 
Charbonnel & Balachandran 2000; Reddy & Lambert 2005; 
Carlberg et al. 2010; Charbonnel & Lagarde 2010; Kumar 
et al. 2011; Martell & Shetrone 2013; Adamów et al. 2014, 
2015; Yan et al. 2018; Li et al. 2018; Deepak & Reddy 2019; 
Gao et al. 2019; Singh et al. 2019). Moreover, ~ 6% of these 
"Li-rich giants exceed the meteoritic abundance of 3.3 dex, 
indicating that additional "Li must have been generated or 
deposited within them (Balachandran et al. 2000; Zhou et al. 
2019; Singh et al. 2019). The engulfment of SBs is a possible 
explanation for high surface "Li abundances (Sandquist et al. 
1998; Siess & Livio 1999; Sandquist et al. 2002; Aguilera- 
Gómez et al. 2016a,b; Soares-Furtado et al. 2021) because 
SBs do not reach the requisite temperatures to burn their pri- 
mordial "Li. However, there are other pathways for lithium 
enrichment, such as the Cameron & Fowler (1971) mecha- 
nism, which acts after the early red giant branch (RGB). The 
existence of these different pathways makes it harder to iden- 
tify the source of enrichment for stars after the early RGB. 
Infrared excess is a potential indicator of stellar mass loss 
from engulfment, and evolved stars with infrared excess tend 
to have increased "Li and rotation rates (Mallick et al. 2022). 

Several analytical studies (Metzger et al. 2012; Yamazaki 
et al. 2017; Jia & Spruit 2018) have focused on planetary en- 
gulfment by main-sequence (MS) or pre-main-sequence stars, 
where envelope ejection is unlikely because of the high gravi- 
tational binding energy. As for post-MS planetary engulfment, 
early analytical estimates (Nelemans & Tauris 1998; Livio & 
Soker 1984) suggest that SBs with masses! S 10.Mjup cannot 
unbind the stellar envelope. Staff et al. (2016) simulated the 
engulfment of a massive planet by stars in the RGB and AGB 
using 3D hydrodynamics. However, their results regarding en- 
velope ejection were limited by numerical resolution. Overall, 
planetary engulfment remains a relatively unexplored problem 
in the context of hydrodynamical simulations. 


l We use the International Astronomical Union nominal values for solar 


system constants (PrSa et al. 2016). 
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Figure 1. Distribution of exoplanets more massive than Jupiter 
(NASA Exoplanet Archive 2022) over the ratio of their geometrical 
and gravitational radii at the onset of engulfment, assuming they are 
engulfed at their current separations. Since tidal decay will lead to 
the engulfment of planets at smaller orbital separations than their 
current ones, this value is likely a lower limit. The geometrical radii 
are larger for 82% of these exoplanets. 


Previous work on CEE has used the “wind tunnel” numeri- 
cal formalism to study the flow in the vicinity of the engulfed 
companion, accounting for the density gradient in the stellar 
envelope (MacLeod & Ramirez-Ruiz 2015a,b; MacLeod et al. 
2017; Murguia-Berthier et al. 2017; De et al. 2020; Everson 
et al. 2020). These density gradients change the flow mor- 
phology and give angular momentum to the gravitationally 
focused gas, thereby changing the drag forces on the compan- 
ion. Most of this previous work has focused on interactions 
between an evolved star and a compact companion, for which 
gravitational drag dominates. For substellar companions, ram 
pressure drag might dominate, depending on stellar structure 
and on the companion. While some studies have recognized 
the importance of ram pressure (Staff et al. 2016; Jia & Spruit 
2018), it has not yet been accounted for in detail. 

The drag forces acting on the engulfed SB influence the 
dynamics of the orbital decay, the observational signatures 
associated with it, and ultimately whether the SB can survive 
engulfment. Here we study planetary engulfment using the 
wind tunnel numerical formalism. 

In Section 2, we discuss the physical processes relevant to 
engulfment, particularly the relative importance of ram pres- 
sure and gravitational drag forces. In Section 3, we provide 
a brief review of the wind tunnel framework, and discuss its 
applicability and limitations in the context of planetary en- 
gulfment. We discuss flow morphology and drag coefficients 
in Section 4, and apply these results to planetary engulfment 
in Section 5. Section 6 summarizes our results. 


2. PHYSICAL ASPECTS OF ENGULFMENT 


2.1. Gravitational and geometrical regimes 
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Figure 2. Dimensionless flow parameters (left: ratio between geometrical and gravitational radii, right: number of density scale heights across 
the larger radii) as a function of substellar body (SB) mass and position inside a model of the Sun evolved to 10.25 (top) and to the tip of the red 
giant branch (bottom). Dashed lines show the transition between geometrical and gravitational regimes, and between mild and strong density 
gradients. Dotted lines show estimates for SB disruption, either by tidal forces (when the separation equals the tidal radius R+) or by ram pressure 


(f — 1, equation 4). 


The interactions of an engulfed SB with the surrounding 
stellar material are gravitational (gravitational drag) and ge- 
ometrical (ram pressure drag on the geometrical surface of 
the SB). Gravitational drag arises from gravitational focus- 
ing of material behind the SB; this focused material exerts a 
force against the direction of motion. An engulfed SB with 
mass Msg travels relative to the surrounding gas at an orbital 
speed Vorb. Gas with an impact parameter smaller than the 
gravitational radius of the SB, 


(1) 


will be gravitationally focused behind the SB. This gas will 
exert a force F} = 1C; R2 pv2, (Hoyle & Lyttleton 1939; 
Bondi & Hoyle 1944; Bondi 1952), where G is the gravita- 
tional constant (we use the value given in Tiesinga et al. 2021), 
p is the envelope mass density, and C, is a dimensionless co- 
efficient of order unity. On the other hand, the pressure field 


Ra = 2GMsg/vÀ,, 


at the surface of the SB will exert a ram pressure force of 
the form Fp = 1C, R25pv2,, where Rsg is the geometrical 
radius of the SB, and C, is a dimensionless coefficient of 
order unity. 

Drag forces are approximately equal to the momentum per 
unit time passing through the cross-section for the correspond- 
ing interaction (geometrical or gravitational). The ratio be- 
tween ram pressure and gravitational forces is therefore equal 
to the ratio of the cross-sections, (Rsg / Ra)”, or equivalently 
the ratio (Vorb [es SB)", where vecsp = 2G Mss / Hag is 
the escape velocity from the SB. 

CEE studies have dealt almost exclusively with the engulf- 
ment of a compact object, such as a neutron star or black hole, 
in which case the interactions between the companion and the 
stellar material are mostly gravitational. Figure 1 shows the 
ratio between geometrical and gravitational cross-sections at 
the onset of engulfment for the known exoplanets, assuming 
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they are engulfed at their current orbital separations. Planets 
are likely to be engulfed at separations smaller than their cur- 
rent separations as a result of tidal decay. Since the Keplerian 
speed is greater at smaller separations, the gravitational radii 
of planets is likely to be smaller at engulfment than it is today, 
and more planets will be in the geometrical regime at the onset 
of engulfment. Equivalently, planets engulfed at earlier stages 
of stellar evolution are more likely to be in the geometrical 
regime, because they must orbit their host star more closely 
to be engulfed during earlier stages. 

Figure 2 shows the same ratio as Figure 1, but as a function 
of SB mass and position inside a Mọ star evolved to 10 R5 
(top panel) and at the tip of the RGB (bottom panel). We 
computed the properties of this star using the Modules for 
Experiments in Stellar Astrophysics (MESA, Paxton et al. 
2011). As the SB dives deeper into the envelope, its interac- 
tions with the gas become increasingly geometrical because 
the Keplerian speed increases inwards. While it is possible 
for the Keplerian speed to decrease inwards if the enclosed 
mass profile is sufficiently steep, the post-MS envelopes we 
consider here are extended enough that the Keplerian speed 
always increases inwards. Therefore, interactions between 
the SB and the surrounding material become increasingly 
geometrical throughout engulfment. 

Figure 2 also shows that more massive SBs tend to be deeper 
in the gravitational regime. Between 107" Mjyp and 10? Mj, 
SB radius varies only between 0.7 Ryup and 1.2 5,5 (we deter- 
mine the radius of each SB using the mass-radius relations 
from Fortney et al. 2007 and Chabrier et al. 2009). As Msp 
increases, Rsg remains approximately constant, but Ra is pro- 
portional to Msg (equation 1), so Rsp/R, is approximately 
inversely proportional to Msg. 

The right panels in Figure 2 shows the number of density 
scale heights across the SB, 


Ep = max (Rsp, Ra) / H,, Q) 


where the density scale height H, = —p/(dp/dr). The 
dimensionless density gradient £, quantifies the heterogeneity 
of the flow at the scale of the SB. Engulfed SBs typically 
experience mild (0 < €, < 1) density gradients. 


2.2. Destruction of the substellar body 
The SB will be tidally disrupted when 


(Penc) X (psp); (3) 


where (penc) is the average density of the material enclosed by 
the SB's orbit, and (psg) is the average density of the SB. This 
criterion is equivalent to the orbital separation being equal to 
the tidal radius of the SB. 

We estimate the SB will be disrupted by ram pressure when 
the kinetic energy per unit volume of the incoming flow equals 
the binding energy per unit volume of the SB, i.e., when (Jia 
& Spruit 2018) 


puža 
(pss) Ue 


=1, (4) 
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Figure 3. Comparison between the binding energy of exterior ma- 
terial Ex, and orbital energy A E>, deposited by substellar bodies 
(SBs) of different masses, as a function of position inside a model of 
the sun at the tip of the red giant branch. Orbital energy deposition 
lines turn dashed at the estimated destruction point for each SB. 


As shown in Figure 2, tidal disruption is the dominant destruc- 
tion process for most SBs. Ram pressure will disrupt SBs 
with masses < 1Mjyp engulfed early in the red giant branch. 


2.3. Envelope ejection 


We use the standard energy formalism of CEE (van den 
Heuvel 1976; Webbink 1984) to estimate analytically whether 
envelope ejection is possible. In this formalism, an engulfed 
companion will eject material exterior to the orbital separation 
r if the binding energy of that material is smaller in magnitude 


than the change in orbital energy of the SB, i.e., |AE*,| > 
| Es], where 
" G Menc Ms GMenc,oM: 
= o 06 
r 2ro 
Mene,0 GM! 
Evina = f (-5— + u) dMine- (6) 
Mene T 


Here, u is the specific internal energy, Menc is the enclosed 
mass at orbital separation r, and the 0 subscript refers to 
values at the initial time t = 0. These two equations neglect 
the binding energy between the star and the SB; at the low 
mass ratios we study here, the envelope is much more bound 
to the core than to the SB, so we omitted these terms for 
conciseness and readability. Figure 3 compares Er. for a 
MESA model of the Sun at the tip of the RGB to AE», for 
SBs of several masses. We define the core-envelope boundary 
as the radial coordinate at which the hydrogen mass fraction 
is 1/10. We consider envelope ejection to be possible if the 
[A E54] is at any point (before destruction of the SB) above the 
value of the | £;.,4| curve at the core-envelope boundary. The 
Figure shows this value as a horizontal line labeled "envelope 
binding energy.” This prescription assumes that energy can be 
efficiently distributed in the envelope, so that energy deposited 
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Figure 4. Ratio of geometrical and gravitational cross-sections at 
the point of envelope ejection or SB destruction, as a function of 
stellar radius and SB mass. The dashed line shows the minimum 
mass required to eject the envelope according to the standard energy 
formalism. When engulfment ends, either by destruction of the SB 
or envelope ejection, all SBs are in the geometrical regime. 


at r can eject material at r’ < r. As we will see (Figure 6), the 
SB significantly disturbs material within a few max (Rsg, Ra) 
of its current location. Tidal forces destroy massive SBs at 
orbital separations comparable to their size (see the lines for 
10Mjyp and 80Mj,, in Figure 3), so we expect at least the 
energy they deposit near the destruction point to reach the 
core-envelope boundary. 

According to the standard energy formalism, objects with 
masses = 10Mjyp might eject the envelope of a Sun-like star 
at the tip of the red giant branch. Figure 4 shows Rsp/Ra 
at the location of either SB destruction or envelope ejection, 
as a function of SB mass and stellar radius. The solid line 


shows the minimum SB mass required to eject the envelope. 


The star expands throughout the RGB, so its binding energy 
decreases and envelope ejection by smaller SBs becomes 
possible. No companion with mass < 80Mjup can eject the 
envelope of an Mo star early in the RGB. All SBs are in 
the geometrical regime when they are destroyed or when they 
eject the envelope, highlighting the need for numerical models 
of ram pressure drag. 


2.4. Orbital decay timescales 


The rate of orbital energy dissipation is the work per unit 
time done by the drag forces, so the orbital decay timescale 
for ram pressure drag is 


Eor G Msp Menc G Menc Msp 
Tinsp,p ~ = = = 2.3 . (7) 
Eob 2a FpVorb 2n Cya Rg V bP 
and for gravitational drag 
GMsp Menc Menc Vorb 
Tinsp.g RI = . (8) 
2aF Vorb Msg / 81GC,ap 


In the geometrical regime, the orbits of more massive SBs 
decay more slowly because they experience approximately 
the same force, but have larger inertia. As before, radius is 
almost constant in mass between 1 Mj, and 100Myup, so the 
change in the geometrical orbital decay timescale as a result of 
changes in SB radius is negligible. In the gravitational regime, 
however, the orbits of more massive SBs decay faster because 
the gravitational cross-section scales as the square of their 
mass, overcoming the inertial term. Whether the orbital decay 
timescale is equation 7 or equation 8 depends on Rsg /Ra at 
the onset of engulfment. 

Tides might dominate the orbital decay of the companion 
in the outer envelope (Stephan et al. 2020), but we do not 
account for them here. Additionally, the drag orbital decay 
timescales are sensitive to the density of the stellar envelope 
at the onset of engulfment because the object will spend most 
of its time in the outer envelope, where drag forces are smaller. 
The orbital and tidal evolution during the post-main-sequence 
set the initial conditions for engulfment. Hydrodynamical 
simulations (Staff et al. 2016) show that the star will overflow 
its Roche lobe. Our calculations neglect these effects by using 
the unperturbed stellar structure throughout engulfment. 

We expect the orbits of all SBs engulfed near the tip of 
the RGB to decay on a timescale comparable to Equation 8 
because they are in the gravitational regime at the onset of 
engulfment (see Figure 2). In earlier stages of stellar evolu- 
tion, however, less massive SBs can be in the geometrical 
regime at the onset of engulfment, and more massive SBs 
in the gravitational regime. Since the gravitational and ram 
pressure timescales have opposite scaling in mass, the orbits 
of intermediate-mass SBs will decay the slowest in these stars 
(see Appendix A). 


3. WIND TUNNEL NUMERICAL FRAMEWORK 


Global simulations of planetary engulfment that account for 
changes in the internal structure of the post-MS star and the 
SB are computationally challenging and expensive because 
they involve two vastly different scales. The scale of the flow 
in the vicinity of the SB is ~ max (Rss, Ra) ~ Rp, while 
the scale of the orbit and the star is ~ 1 au ~ 10? Ry. This 
disparity of scales motivates isolating the processes occurring 
at each scale, not only for computational feasibility and accu- 
racy, but also to understand the role of each of these processes 
and eventually the interplay between the processes at different 
scales. Here we perform simulations of the flow within a few 
max (Rsg, Ra) of the SB. We aim to understand the mor- 
phology of the flow in the vicinity of the SB, and the forces 
acting on it. Measurements of these forces allow numerical 
integration of the equation of motion of the SB inside the star. 
This approach allows inexpensive yet reasonably accurate 
exploration of a much larger region of parameter space. 

We use the “wind tunnel” numerical setup (MacLeod et al. 
2017), illustrated in Figure 5. We simulate a local domain 
in the frame of the engulfed SB, and we supply a time- 
independent “wind” from the —ĉ direction. The flow mor- 
phology (and therefore the drag coefficients) is uniquely de- 
termined by a set of dimensionless parameters: the mass 
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ratio q = Msp/Menc, the dimensionless density gradient £, 
(equation 2), the mach number M, and the ratio between 
geometrical and gravitational radii Rsg/Ra. The equation 
of hydrostatic equilibrium for the envelope (MacLeod et al. 
2017) relates these parameters as 


. Ra EN 2q —4 2 
min (i a) Ep = xg = ay fe M, (9) 


where fp is the ratio between the speed of the companion 
relative to the surrounding material and the Keplerian speed. 
Equation 9 implies three dimensionless parameters are enough 
to uniquely specify the flow morphology; for most of this work 
we use the set (q, €p, Rsp/Ra). Since the drag coefficients are 
a function of only these dimensionless parameters, a single 
simulation can determine the drag for a variety of physical 
systems. The drag coefficients do not depend on dimensional 
quantities that set the scale of the system. These quantities 
do not change the flow morphology, and the total drag scales 
with them in known ways (for example, linearly, in the case 
of the density). For simplicity, we set the density at y — 0 to 
1gcm-?, and the wind speed to 1 cm s^ +. We set the pressure 
at y = 0 so that the flow has the Mach number implied 
by hydrostatic equilibrium and the rest of the dimensionless 
parameters. 

After computing the flow properties at y — 0, we integrate 
the equations of hydrostatic equilibrium for a massless at- 
mosphere up to the boundaries in the +y directions. In the 
—y direction, we extend hydrostatic equilibrium to the ghost 
zones. We add the gravitational field of the enclosed stellar 
mass so that the gas remains in hydrostatic equilibrium in the 
absence of external forces. The external force in our simu- 
lations, leading to the deflection of the gas, is the gravity of 
the SB. For all other boundaries, we use outflow boundary 
conditions, in which gas can leave the domain but not enter it. 
For more details, see MacLeod et al. (2017). The numerical 
setup is publicly available (see Appendix C). 

We set fk = 1, thereby assuming a circular orbit and no 
corotation between the star and the SB. During the post- 
main-sequence, synchronization tides shrink the orbit of the 
companion and transfer its angular momentum into the star. 
Small SBs do not have enough orbital angular momentum 
to bring the star into corotation (equality of the orbital and 
rotational periods of the star), leading to orbital decay through 
the Darwin (1879) instability. The ratio between rotational 
and orbital periods at the onset of the engulfment of a Jupiter- 
like planet by a Sun-like host star evolved to 0.05Rọ is ~ 
0.2 (Gallet et al. 2017); for gas giant planets engulfed by a 
1.5Mo host star, this ratio is < 0.3 (Table A.1 of Privitera 
et al. 2016b). More massive SBs, such as brown dwarfs, will 
enhance stellar rotation more. 

In the wind tunnel framework, the density gradient (point- 
ing in the — 4 direction) and the velocity of the stellar material 
in the frame of the SB (pointing in the ĉ direction) are perpen- 
dicular. This configuration doesn't accurately represent the 
flow around the SB when its orbit is eccentric. Tides signifi- 
cantly lower the orbital eccentricity of closely orbiting planets 


Figure 5. Schematic of the wind tunnel numerical setup. Left: 
substellar body (SB) embedded in the envelope of a giant star, with 
an orbital decay trajectory shown as a dashed line. Right: density 
slice and velocity streamlines in a wind tunnel simulation of the 
embedded object. The wind velocity and density gradient point in 
the and — 4 directions, respectively. The orbital speed is higher 
than the local sound speed, so the object will drive a shock as it 
moves through the stellar envelope. 


throughout stellar evolution, while the eccentricity of distant 
planets remains roughly constant (see Figure 8 of Villaver 
et al. 2014). The average eccentricity of the known exoplan- 
ets more massive than Jupiter around stars between 0.7 M. 
and 1.5 Mg is ~ 0.2 (NASA Exoplanet Archive 2022). Asa 
result of tides, these planets will likely be engulfed in orbits 
more circular than their current ones. However, there is theo- 
retical and observational evidence for a transient population 
of exoplanets orbiting evolved stars at moderate eccentric- 
ity (Villaver et al. 2014; Grunblatt et al. 2022), whose orbits 
decay before they circularize. 

This framework also neglects the curvature of the veloc- 
ity field within the domain. This approximation is valid 
when the half-length of the domain Ldomain/2 is much smaller 
than the orbital separation r. We set the domain length to 
10 max (Rsg, Ra). In the geometrical regime, the condition 
Laomain/2 < r reduces to Rsp/Ra < (1 + q) /10. In the 
gravitational regime, it reduces to q « 1/9. Figure 8 shows 
these constraints. 

We wrote this numerical framework as a setup for the 
FLASH (Fryxell et al. 2000; Dubey et al. 2014, 2015) code. It 
uses FLASH to solve the equations of inviscid hydrodynamics 
on a Cartesian mesh with adaptive mesh refinement. We use 
an ideal gas equation of state P = (y — 1) u, where P is the 
pressure, ~ the ratio of specific heats (which we take to be 
5/3), and u is the internal energy per unit volume. The base 
resolution of our simulations is 160 cells per dimension. We 
use adaptive mesh refinement with a criterion based on prox- 
imity to the surface of the SB; we choose the maximum level 
of refinement such that there are always ~ 60 cells across the 
SB. See Appendix A.1 for hydrodynamics convergence tests. 


3.1. Model for the substellar body 


Previous wind tunnel CEE simulations between extended 
stars and compact objects have modeled the compact object as 
a numerical “sink.” Inside the sink, the simulation multiplies 
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Figure 6. Density (in units of upstream density at y = 0, po.) slices in wind tunnel simulations, for several values of the ratio between geometrical 


and gravitational radii Rss / Ra, and the number of density scale heights per max (Rsg, Ra). These slices are at t = 20 max (Rsg, Ra) / vow, 


where Vors is the upstream gas speed. At low Rsp/ Ra, a spherically symmetric envelope of material that cannot be accreted forms around the SB, 


suppressing ram pressure drag. At high Rss / Ra, a vacuum forms behind the SB, increasing ram pressure drag and suppressing gravitational 
drag. An animated version of this figure is available in the HTML version of the article. The animation shows the time evolution of the density 


slices from t = 0 to t = 20 max (Rsg, Ra) / vo. 


fluid variables by a small number, creating a numerical vac- 
uum that emulates accretion of the surrounding material onto 
the object. In real systems, these objects accrete because the 
material accumulating around them is hot and dense enough 
to cool via neutrino emission (MacLeod & Ramirez-Ruiz 
2015b; Fragos et al. 2019). There is no such cooling channel 
for material near the surface of the SB, and the gas is too 
optically thick to cool radiatively. The timescale over which 
radiation will carry energy through the optically thick sur- 
rounding material, allowing it to be accreted, is much longer 


than the orbital decay timescale. We therefore model the SB 
as a rigid body with a reflective surface. We use FLASH’s 
unsplit hydrodynamics solver, which is based on the unsplit 
magnetohydrodynamics solver (Lee 2013). The code applies 
the reflective boundary condition at the surface of the rigid 
body, i.e., in the rigid body cells adjacent to the fluid cells. 
The reflective boundary implementation in FLASH requires a 
Courant-Friedrichs-Lewy (CFL) number < 0.3 for numerical 
stability because the reconstruction order near the boundary 
is lower. 
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Figure 7. Ram pressure (left) and gravitational (right) drag coefficients for a set of simulations with mass ratio q = 2.15 x 107°, as a function 


of the ratio between geometrical and gravitational radii, Rss /Ra. Each line corresponds to a different dimensionless gradient €p. The drag 


coefficients depend on Rsg / Ra most strongly during the transition regime Rss ~ Ra, then have a weaker dependence until the mach number 


approaches (and goes below) unity. 


During the engulfment of a compact object by a giant star, 
the geometrical size of the object is negligible (Rsg < Ra), 
whereas an engulfed SB can have Rsg = Ra. When the 
radius of the “sink” object is = 0.2R,, the shock morphol- 
ogy changes into a "tail shock" that trails the path of the 
object (Figure 10 in Ruffert & Arnett 1994). As we will see 
(Figure 6), the reflective boundary prevents this change in 
morphology. These qualitative differences motivate modi- 
fied wind tunnel simulations that more accurately represent 
how engulfed SBs interact with their surroundings, and the 
associated flow morphology and drag coefficients. 

We do not model the internal structure of the SB, whose 
mass loss and deformation could affect the flow morphology 
around it, and the cross-section for interactions with its sur- 
roundings. Ram pressure will flatten the surface of the SB 
facing the incoming flow and compress the SB, making it 
harder to disrupt (Jia & Spruit 2018). Gradual ablation of the 
SB is unlikely to destroy the SB before “global” processes 
that act on its dynamical timescale (Passy et al. 2012; Jia & 
Spruit 2018), although the hydrodynamics of ablation in this 
context are uncertain. We discussed these global processes 
in Section 2.2; we evaluate their corresponding destruction 
criteria using only the unperturbed SB structure. 


3.2. Drag force measurements 


We measure the forces on the object when the simulation 
reaches steady state, which takes a few flow-crossing times 
7, = max (Ra, Rsg) /Vorn. We measure the gravitational drag 
force by integrating the gravitational force of the surrounding 
density field up until 1.6 max (Rsg, Ra), as in MacLeod et al. 
(2017). The ram pressure drag force is 


F, = -4 PdA, (10) 
S 


where dA is the area element and S is the surface of the 
SB. Our setup uses FLASH to solve the equations of invis- 
cid hydrodynamics; while the discretization of the equations 
results in numerical viscosity, there is not a boundary layer 
around the SB, which could change the ram pressure drag 
and the viscous stresses acting on it. For the same reason we 
do not study the dependence of the ram pressure drag on the 
Reynolds number. 

The coefficients measured from the steady state simulations 
are valid if the timescale over which the simulation reaches 
steady state is much shorter than the orbital decay timescale. 
For SBs dominated by ram pressure, whose orbital decay time 
is equation 7, Tc/Tinsp.p = p/(psB) < (penc) /(Psp), Where the 
inequality holds because the density decreases monotonically 
with radius. This ratio is less than unity at the onset of en- 
gulfment; when it reaches unity, the core tidally disrupts the 
SB (equation 3). For SBs dominated by gravitational drag, 


Tel Tinspg = (Msp / Max)? D/ (penc) <1. 


4. FLOW MORPHOLOGY AND DRAG 
4.1. The gravitational and geometrical regimes 


Figure 6 shows steady state density slices of wind tunnel 
simulations with q = 2.15 x 1072, at several Rsg/Ra and Ep 
At low Rsg /Ra (the leftmost column in the Figure), the SB 
gravitationally focuses gas as in the Bondi-Hoyle-Lyttleton 
formalism. However, since the SB cannot accrete, this focused 
material accumulates at its surface. The pressure force exerted 
by this material opposes the pressure force from material ac- 
cumulated in front of the object a result of compression when 
the SB moves through the gas. The resulting pressure field at 
the surface of the object is spherically symmetric, suppressing 
ram pressure drag (Thun et al. 2016). As Rsp/Ra decreases, 
ram pressure drag becomes less important not only because 
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Figure 8. Trajectories in parameter space followed by substellar 
bodies (SBs) of varying masses after being engulfed by a 1Mo 
star at varying evolutionary stages. Each line style represents a 
different evolutionary stage (see legend on top panel), and each color 
represents an SB of a different mass (see inline labels on bottom 
panel). Each line starts at a radial coordinate r = 0.9R,, and ends 
at the point of SB destruction, represented by a triangle. The squares 
represent the hydrodynamical simulations done in this work. 


the geometrical cross-section decreases, but also because of 
morphological changes to the flow. 

When Rsg ~ Ra, the gravity of the SB is not strong enough 
to deflect the surrounding material, and a low-density region 
forms behind the SB. The ram pressure force exerted by gas 
in front of the SB is now unopposed, and gravitational drag 
decreases significantly. When Rsg >> Ra, the material in 
front of the SB dominates gravitational drag (as opposed to 
material behind the SB, as when Rsg < Ra), and the gravita- 
tional drag coefficient becomes negative. Similar results have 
been found for compact objects with outflows (Gruzinov et al. 
2020; Li et al. 2020; Kaaz et al. 2022) and luminous planetesi- 
mals moving through a disk (Masset & Velasco Romero 2017; 
Masset 2017). In those settings, feedback from the object 
can interfere with the flow at impact parameters ~ Ra that 
the SB would have gravitationally focused had there been no 
feedback. In our simulations, the rigidity of the SB and its 
reflective surface have this effect when Rsg > Ra. 

Figure 7 shows the drag coefficients for the same set of 
simulation parameters as Figure 6. In the gravitational regime, 
the drag force is most naturally written as 


R 2 
Fas, = -T pv? R? c,+ (=) C,|9, (D 


where 6 is the unit vector in the direction of the SB's velocity. 
On the other hand, in the geometrical regime 


R —2 
Fas; = —Tpv" Rp c+ RE) C lô aD 


Equations 11 and 12 are equivalent, but it is often convenient 
to write the drag force in the functional form of the dominant 
source of drag, with a small correction term that accounts for 
the nondominant source of drag. The drag can be written as a 
combination of the previous two expressions, 


Farag = —rpv? max (Rss, Ra)” (Cre + Coe) 0, (13) 


where 
IRE 
Cp,ett = Cy min | 1, EJ i (14) 
, Ra 
Rl 
Cg eft = Cg min | 1, | | (15) 
Rsg 


The primary motivation for writing the drag force this way is 
that the ratio between the “effective” coefficients equals the 
ratio between the drag forces. Additionally, when the flow con- 
verges to the geometrical regime, the effective gravitational 
drag coefficient approaches zero. The gravitational drag coef- 
ficient alone does not approach zero when Rsg > Ra because 
in that limit the material in front of the SB dominates gravita- 
tional drag. Since larger SBs have more material in front of 
them, the gravitational force exerted by that material is larger. 
However, this gravitational drag force increases with Rsg /Ra 
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Figure 9. Density (in units of upstream density at y = 0, po.) in wind tunnel simulations as a function of space. Each panel shows a different 
combination of the dimensionless density gradient (€p), and the ratio between the mass of the substellar object and the mass enclosed by its orbit 
(q). The ratio between geometrical and gravitational radii is 0.3 in all panels. These slices are att = 20Ra / Vorb, Where Vorb is the upstream gas 
speed. Larger density gradients and Mach numbers result in higher drag coefficients. An animated version of this figure is available in the HTML 
version of the article. The animation shows the time evolution of the density slices from t = 0 to t = 20Ra/Von. 


slower than quadratically, so gravitational drag becomes a 
smaller fraction of the total drag as Rsg /Ra increases. 

Figure 7 shows that when Rss < 0.25 the effective drag 
coefficients are independent of Rsg /Ra, so that 


Cg eff == Cg (9, €p) , (16) 
Cy ett = Cp (Rsp/ Ra)” ex 0. (17) 


On the other hand, the ram pressure drag coefficient de- 
pends on Rsp/R, even when Rsg >> Ra. This dependence 
can be understood from the equation of hydrostatic equilib- 


rium, which when Rsg > Ra reduces to 


Ra 2q 


= 3 Mf. 18 
Ra” Gag? pm 


At a fixed mass ratio q and dimensionless density gradient 
Ep, the mach number M œ (Rss/Ra) JT. As Rsp/Ra 
increases, the mach number decreases, so the density disconti- 
nuity across the shock is smaller, reducing ram pressure drag. 
This gradual transition towards subsonic flow can be seen in 
the top row of Figure 6, from left to right. In the geometrical 
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Figure 10. Orbital decay trajectory of a 1 Myup planet inside a 1Mo 
star evolved to 10Ro, as computed using analytical (left panel) and 
numerical (right panel) drag coefficients. Trajectories turn transpar- 
ent at the estimated point of destruction. Orbital decay takes ~ 24d 
(a factor of ~ 2 faster) when using numerical drag coefficients. 


regime, 


C, et = Cy (Ra/Rsp)” z 0, (19) 
Cuan = Cp = Cp (q, €p, Raps Ra). (20) 


4.2. Dependence on mass ratio and dimensionless density 
gradient 


The dependence of the drag coefficients on the density gra- 
dient, the Mach number, and the mass ratio can be understood 
from the relationship between these parameters in hydrostatic 
equilibrium (Equation 9), and the flow morphology in Fig- 
ure 9. Simulations with stronger density gradients have larger 
drag coefficients because the SB interacts with higher den- 
sity gas (from deeper in the envelope) both geometrically and 
gravitationally. From hydrostatic equilibrium, increasing the 
mass ratio at a fixed Mach number will result in a stronger 
density gradient, and larger drag coefficients. 

Similarly, at a given density gradient, larger Mach numbers 
result in a narrower shock opening angle, allowing focused 
material to accumulate closer to y = 0, increasing the hori- 
zontal component of the drag force and therefore the gravi- 
tational drag coefficient. Higher Mach numbers also result 
in increased gas compression in front of the SB, increasing 
the ram pressure drag coefficient. At fixed density gradient, 
increasing the Mach number requires decreasing the mass 
ratio, so increasing mass ratios reduce the drag coefficient. 


5. APPLICATIONS TO PLANETARY ENGULFMENT 
5.1. Parameter space for hydrodynamical simulations 


Figure 8 shows the trajectories of engulfed SBs in parameter 
space, for several SB masses and stellar evolutionary stages. 
Each line style corresponds to a stellar evolutionary stage, 
and each color corresponds to an SB mass. Each line starts 
at a radial coordinate r = 0.9R, and ends at the point of SB 
destruction, represented by a triangle. This figure, as Figure 2, 
assumes SBs are in circular orbits throughout engulfment. 
The wind tunnel framework cannot simulate the region in 
parameter space where the domain size would be compara- 
ble to the orbital separation (the gray region in the top and 
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Figure 11. Orbital decay trajectory of a 80Mjup brown dwarf inside 
a 1 Mg star evolved to 100Ro, as computed using analytical (left 
panel) and numerical (right panel) drag coefficients. For readability, 
we show the trajectories only down to an orbital separation of 2.5Ro. 
While the orbital decay timescale is similar in both cases (z 1.3 yr), 
numerical drag coefficients result in a more eccentric orbit. 


middle panels). We discussed this limitation quantitatively in 
Section 3. 

The top panel shows that more massive SBs (those with 
higher q) are deeper in the gravitational regime, and that 
Rsp/Ra increases throughout engulfment. The destruction 
points for SBs with q = 5 x 10^? lie in a line; these SBs are 
destroyed by tidal disruption. From the definition of the tidal 
radius, tidal disruption will occur when the orbital separation 
r = Rggq- "3. Under the assumption of circular orbits, Ra = 
2qr. Combining these two equations, Rsg/ R, = q- ?/? /2 at 
the point of destruction, as seen in the figure. 

This figure determines the parameter space that we must 
simulate to capture the diverse flow morphologies of the sys- 
tems that undergo engulfment. In Section 4.1, we found that 
the drag coefficients are independent of Rsg /Ra in the gravi- 
tational regime. The blue region in Figure 8 shows this regime. 
Therefore, hydrodynamical simulations need to span only the 
transition and geometrical regimes. Each square in the Figure 
represents a hydrodynamical simulation, of which we ran 428. 


5.2. Orbital decay trajectories 


The equation of motion for the SB inside the star is 


Mss = (Msp — pVes)E- Fam. CD 
where v is the velocity of the SB, t is time, Vsg is the vol- 
ume of the SB, and g = —GMenci/ r? is the gravitational 
acceleration from the mass enclosed by the SB's orbit. The 
term pVsp& is the buoyancy acting on the SB, which is im- 
portant when the local density equals the average density of 
the SB. Since the local density is always smaller than the 
average enclosed density, and the SB will be tidally disrupted 
approximately when the average enclosed density is equal to 
its own mean density, buoyancy won't be important before 
SB destruction. 

We integrate equation 21 numerically using the IAS15 (Rein 
& Spiegel 2015) integrator from the N-body code REBOUND 
(Rein & Liu 2012). At every timestep, we compute the dimen- 
sionless parameters (q, €p, Fisp / Ra) and linearly interpolate 
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Figure 12. Top left panel: orbital decay and energy transport timescales as a function of the radial coordinate r of the substellar body (SB) inside 
the star. The orbital decay timescales for the two more massive substellar bodies show peaks because we define the orbital decay timescale 


as |r/7|, and as a result of eccentricity in the orbit 7 changes sign. Bottom left panel: rate of energy deposition and luminosity (excluding the 


stellar luminosity), again as a function of position inside the profile. When the energy transport timescale time is much longer than the orbital 


decay timescale, the luminosity is much smaller than the energy deposition rate. Right panel: luminosity as a function of time, in units of the 


unperturbed stellar luminosity. 


the drag coefficients C; and Cp we measured in the hydro- 
dynamical simulations. For points outside the domain of the 
interpolation, we used the nearest available point. 

We initially place the SB in a circular orbit at an orbital 
separation 0.9R,. At every timestep we interpolate the prop- 
erties of the stellar profile using the GSL implementation 
of the Steffen method (Steffen 1990). We stop the integra- 
tion when r < Rsg, and apply the destruction criteria during 
postprocessing. 

We compute stellar profiles using MESA. The orbital decay 
timescale is < 400 yr for the systems we consider here (a 
1 Mjup planet inside a 1 Mo star at the tip of the RGB takes ~ 


400 yr). On the other hand, the timescales R/R and L/L over 
which radius and luminosity change significantly, respectively, 
are > 10° yr throughout the RGB. Therefore, stellar structure 
doesn't change significantly as a result of stellar evolution 
over the timescales relevant to engulfment. We do not model 
the effects of engulfment on stellar structure, and use a single 
unperturbed MESA profile throughout our integration of the 
equation of motion of the SB. 

Figure 10 and Figure 11 contrast the orbital decay trajecto- 
ries obtained by using either analytical (Cj, = 0.25, C, = 1) 
or numerical drag coefficients, for two different systems. We 
chose systems whose orbital decay timescales were the small- 
est compared to the orbital period at their initial separations, 
so that their trajectories were easy to visualize. Figure 10 
shows the trajectories for a Jupiter-like planet inside a 1M, 
star evolved to 10 Ro. The orbit remains nearly circular, so 


the energy deposition profiles are similar when using ana- 
lytical or numerical drag coefficients. However, the orbit 
decays a factor ~ 2 faster with the numerical coefficients. 
The timescale of energy transfer from the orbit into the star 
determines whether the star will transport the energy to the 
surface, or react dynamically. This distinction is particularly 
relevant for the stellar envelope, in which convection can 
carry energy to the surface quickly, lowering the efficiency of 
energy deposition (Wilson & Nordhaus 2019, 2020, 2022). 

Figure 11 shows the trajectories for a 800Mj,5 brown dwarf 
inside a star evolved to 100 R5. For this system, both models 
for the drag coefficients yield similar orbital decay timescales, 
but the orbit is significantly more eccentric when using nu- 
merical drag coefficients. As discussed in Section 3, our 
hydrodynamical simulations neglect eccentricity, and there- 
fore do not capture the flow morphology when the orbit of 
the SB is significantly eccentric. However, they show that at 
least small eccentricities arise from initially circular orbits 
faster than suggested by analytical drag. The evolution of 
the eccentricity is important for many types of systems. For 
the most massive SBs that can eject the envelope, it deter- 
mines the eccentricity of their orbit around the stellar remnant 
(Szólgyén et al. 2022). For small SBs, the eccentricity at 
the point of destruction determines the redistribution of their 
enriched material throughout the star. 


5.3. Stellar brightening 
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Figure 13. Engulfment outcomes for known substellar bodies (SBs, 
NASA Exoplanet Archive 2022), as a function of SB mass and or- 
bital separation. This figure assumes that SBs are engulfed at their 
current separations and that they all orbit 1Mo stars. At a fixed SB 
mass, envelope ejection becomes easier as the star ascends the red 
giant branch because its binding energy decreases. A large fraction 
of SBs engulfed early in post-main-sequence evolution will dissolve 
in the convective region, yielding potential surface abundance en- 
hancements of the "Li isotope. We also show systems with orbital 
separations larger than the maximum stellar radius during the RGB; 
if tidal forces lead to the engulfment of massive brown dwarfs with 
orbital separations larger than this maximum RGB radius (1.e., the 
brown dwarfs in the upper right region of the figure), they could eject 
the stellar envelope. 


Planetary engulfment results in a transient dominated by 
recombination in mass ejected from the outer layers of the 
star (Metzger et al. 2012). Since we do not model changes 
to stellar structure, we study only the long-term emission 
resulting from the eventual transport of orbital energy from 
the deep layers of the star to its surface. If an energy dE is 
deposited into the star, the upper bound (assuming all energy 
leaves as radiation) for the average increase in stellar lumi- 
nosity is dL = dE /Tgu, where ky is the energy transport 
timescale from the location of energy deposition to the sur- 
face. More generally, for continuous energy deposition, the 
time-averaged additional luminosity at time t is 


- 
AL (t) = I JE d (22) 
to 


TKH 


where we determine tọ by noticing that energy deposited 
a time t — to ago contributes to the increase in the average 
luminosity only if t— to < Tpu (to), since if t—to > Tpu (to), 
the energy deposited at to has already been radiated away. If 
the energy deposition increases sharply in the inner regions 
of the star, E ~ 6 (t/ — t) AE (t^) in equation 22, and AL ~ 
A E /txu (e.g., MacLeod et al. 2018). 


We compute E from the work per unit time done by the 
drag forces, 
E = —Firag ' V. (23) 


We compute the energy transport timescale by adding the cell 
crossing times in the stellar profile. The most effective energy 
transport mechanism at each cell (either radiative diffusion or 
convection) determines the crossing time for that cell. 

Figure 12 shows the quantities involved in our calculation 
of the average luminosity for the engulfment of SBs of varying 
masses by a 1 Mo star evolved to 10Rọ. The top left panel 
compares the orbital decay and energy transport timescales as 
a function of the orbital separation (equivalently, time). When 
the orbital decay timescale is much longer than the energy 
transport timescale, the luminosity is close to the energy depo- 
sition rate, whereas deep in the envelope the energy transport 
timescale is much longer than the orbital decay timescale, 
making the luminosity much smaller than the energy depo- 
sition rate. The bottom left panel shows this behavior. The 
right panel shows the luminosity. 

From R, = 10, to the tip of the RGB, the additional 
luminosity from engulfing SBs of most masses 2 My is 
comparable or larger than the stellar luminosity, in some cases 
by several orders of magnitude (see Metzger et al. 2017). SBs 
engulfed early in the post-main-sequence have higher energy 
deposition rates (because of the higher stellar density), and 
the stars that engulfed them are dimmer. On the other hand, 
evolved stars are sparser and more luminous, but have smaller 
energy transport timescales at the point of SB destruction. 
While these luminosity estimates rely on Equation 22 (which 
doesn’t accurately describe radiative transfer inside the star) 
and depend on uncertain processes (such as the destruction of 
the SB), they suggest that engulfment is energetically signifi- 
cant throughout the post-main-sequence. 

After SB destruction, the timescale for the luminosity to 
return to its original value is roughly the energy transport 
timescale at the point of destruction. The energy transport 
timescale at the point of destruction becomes shorter as the 
star ascends the RGB. For a given star, more massive SBs, 
which tend to also be denser as a result of the mass-radius 
relation, will result in a longer increase in luminosity because 
they survive deeper into the envelope. For a Sun-like star 
evolved to 10 Ro, the time it takes for the star to return to its 
original luminosity ranges from ~1 yr for a 1 Mjyp planet to 
~5000 yr for a 80Mjyp brown dwarf, as shown in Figure 12. 
On the other hand, for a model of the Sun at the tip of the 
RGB, the time ranges from ~1 yr to ~25 yr for the same 
range of SB masses. 


5.4. Engulfment outcomes 


Figure 13 shows known SBs (NASA Exoplanet Archive 
2022) as a function of their mass and present-day orbital 
separation. The dashed line near the top-right corner shows 
the minimum SB mass for envelope ejection, assuming all 
deposited energy contributes to the ejection. Since the trajec- 
tories of some of the massive SBs that can eject the envelope 
are likely eccentric (Figure 11), we compute this line using 
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the analytical values for the drag coefficients. For destroyed 
SBs, the figure shows whether they'll be destroyed in the 
convective zone or below it, according to the analytical de- 
struction estimates of Equations 3 and 4. This figure assumes 
SBs are engulfed at their present-day orbital separations, and 
that all SBs orbit 1 Mo stars (the average stellar mass reported 
for these SBs’ planetary systems is 1.12Mọ, with a standard 
deviation of ~ 17%). 

Under these assumptions, Figure 13 suggests massive SBs 
can eject the envelopes of evolved Sun-like stars through 
the transfer of orbital energy. This figure also shows that 
a substantial fraction of known SBs might be destroyed in 
the convective region of their host stars, particularly those at 
closer orbital separations because they are engulfed when the 
star is more compact and disrupts them more easily. The "Li 
contained in these SBs could be carried via convection to the 
surface, resulting in enhanced surface abundances. However, 
the mean molecular weight of the SB's enriched material 
is much higher than that of its surroundings. Some of that 
material could settle in a layer near the base of the convective 
region and eventually reach the radiative core (Vauclair 2004; 
Jia & Spruit 2018). Intermediate-mass SBs (those in the 
central region of the plot) are massive enough to survive 
below the base of the convective region, perhaps resulting in 
opacity changes detectable through asteroseismology. 


6. CONCLUSIONS 


We studied the engulfment of substellar bodies (SBs) by 
evolved stars using hydrodynamical simulations of the flow 
in the vicinity of an engulfed SB (the “wind tunnel” frame- 
work, schematically depicted in Figure 5). The steps in our 
numerical framework are: 


1. Determine the hydrodynamical parameter space for 
planetary engulfment. In particular, the range of val- 
ues for the dimensionless parameters that affect the 
morphology of the flow around an engulfed SB. 


2. Run hydrodynamical simulations that span this parame- 
ter space, characterize the resulting morphologies, and 
measure the drag coefficients for the drag forces acting 
on the SB. 


3. Use the drag coefficients to integrate the equation of 
motion of an engulfed SB, and estimate observational 
signatures and outcomes of engulfment. 


Some of our main findings are: 


* The interactions of engulfed SBs with their environ- 
ment are geometrical (ram pressure drag) and gravi- 
tational (gravitational drag). Geometrical interactions 
become increasingly important throughout engulfment 
(Figure 2, Figure 8). All SBs are in the geometrical 
regime when they are destroyed or when they eject the 
envelope (Figure 4). 


According to semi-analytical estimates, tidal disruption 
is the dominant destruction process for SBs with masses 
= 1Mbup at most stages of stellar evolution. 


* Many SBs will be in the gravitational regime at the 
onset of engulfment, and transition to the geometrical 
regime as their orbit shrinks. Our hydrodynamical sim- 
ulations characterize this transition regime qualitatively 
(Figure 6) and quantitatively (Figure 7). 


The engulfment of an SB could increase the luminosity 
of a star by up to several orders of magnitude. The time 
is takes for the star to return to its original luminosity 
depends on its evolutionary stage and the mass of the 
SB (Section 5.3; Figure 12). 


Massive SBs could eject the stellar envelope via trans- 
fer of orbital energy, if the transferred energy can be 
efficiently distributed within the envelope (Figure 13). 
Small SBs are destroyed in the convective region. 


We discussed the applicability of the wind tunnel frame- 
work to planetary engulfment in Section 3. As implemented 
in this work, the framework assumes that the density gradient 
in the stellar envelope and the velocity of the SB are perpen- 
dicular, and that the SB moves at approximately the circular 
Keplerian speed. These assumptions do not hold for some 
massive SBs, whose orbits develop significant eccentricities 
during engulfment. For this reason, we used analytical drag 
coefficients when computing the minimum mass required for 
envelope ejection. The evolution of the internal structure of 
the SB and the star remains a significant uncertainty in our 
calculations. We model the SB as a rigid body with a re- 
flective boundary; while we approximately account for tidal 
and ram pressure disruption to determine the location in the 
star where the SB will be destroyed, only hydrodynamical 
models can describe these processes in detail. Future work 
could study the evolution of the SB's internal structure to 
determine the conditions and timescales associated with its 
destruction. Simulations of the entire star can help understand 
its response to engulfment. Here we used a simplified model 
for energy transport to estimate the engulfment luminosity; 
more sophisticated stellar models including radiative transfer 
could constrain these observational signatures. The numeri- 
cal framework we introduced here can be used to study the 
dynamics of planetary engulfment using comparatively inex- 
pensive simulations that capture the physics of the flow at the 
scale of the SB. 
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APPENDIX 


A. NUMERICAL TESTS 
A.1. Wind tunnel 


Figure 14 shows the ram pressure drag coefficient in a 
simulation with €, = 1 and Rsg /Ra = 0.6 at several reso- 


lutions. We chose a value of Rsp/R, for which computing 
the ram pressure drag would be the hardest numerically. If 
Rsg > Ra, the SB is a larger fraction of the domain size 
(Laomain = 10 max (Rsg, Ra)), so it is easier to resolve. On 
the other hand, if Rsg < Ra, the pressure field around the SB 
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Time [Ry/Vo0] 
10 


— — 15 cells per Rsg 


Time-averaging 
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Ram pressure drag coefficient 
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— — 123 cells per Rsg 
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Figure 14. Ram pressure drag coefficient as a function of time, 
for a simulation with Rsg/Ra = 0.6 and e, = 1. Each line cor- 
responds to a different resolution at the surface of the object. We 
compute the drag coefficient for a simulation by time-averaging the 
time-dependent drag coefficient between 15 max (Rsg, Ra) /vos 
and 20 max (Rsg, Ra) /vos. 


Orbital decay timescale [yr] 


0.00 10° 10 102 


Companion mass [Mjup] 


Figure 15. Orbital decay time for substellar bodies as a function of 
their mass, inside a model of a 1 Mo star evolved to 10 R5. 


is spherically symmetric (Section 4). The hardest simulations 
in which to measure ram pressure drag are those with inter- 
mediate values of Rsp/Ra, for which Rsg is small compared 
to the length of the domain, but for which the pressure around 
the SB is still asymmetric. The maximum resolution in our 
test has 122 cells per Rsg. The ram pressure drag coefficient 
of a simulation with 31 cells per Rsg had a relative error 
(compared to the simulation with the highest resolution) of 
& 1.596. We set the refinement in all our simulations such 
that there are at least 31 cells per Psp. 
A.2. Equation of motion integration 


We tested the implementation of the two effects we account 
for that deviate from the standard two-body problem: the 


change of the mass enclosed by the orbit of the SB, and the 
drag forces. 

We integrated the orbit of a 1 Mj, companion inside a star 
of radius R, = 10Re without drag forces. The initial radial 
coordinate of the SB is 9Ro, and its initial speed is 70% of 
the circular Keplerian speed. After 10^ dynamical times of 
the star, the specific orbital energy of the companion 


v? GMenc Ms dM! 
Forb,sp zd a +G —— (Al) 
2 T Mene T 


was conserved to within a fractional error 5 x 1077. 

We then integrated trajectories inside the same star for SBs 
of different masses, including drag forces (with drag coeffi- 
cients C, = 1 and C, = 0.25). We initially place the SB 
in a circular Keplerian orbit at 95. For each SB mass, we 
computed the corresponding SB radius using the mass-radius 
relation, as discussed in Section 2. Figure 15 confirms that, 
for this star, in the low and high mass limits the orbital de- 
cay timescale scales as Msg / R2, and Mg, respectively (see 
equations 7 and 8). The work done by drag forces was equal 
to the change in orbital energy to within a fractional error 
2x 10-5. 


B. FITTING FORMUL/E 


We fit the results of Figure 13. The minimum companion 
mass required to eject the envelope of a 1 M, star as a function 
of its radius is 


Msp — 1.35892? + 2.11822 — 196.53 
My — 5.0709 x 10-423 + 0.032124z? — 0.516382 + 1’ 
(B2) 
where x = R,/Rg. This formula is valid when 11.45 < 
R,/Ro € 193.7 and agrees with the Figure to within 1.1%. 
The minimum mass for a companion to survive below the 
base of the convective zone is 


5 
Msg 3 2 —1 j 

— (b b 1 > ix B 
My (baa + bx + bix + ) 20 (B3) 


where x = log), (R,/ Ro), bı = —3.01905, b = 2.89041, 
b4 = —0.767818, and 


ai = (5.7113, —23.548, 36.3736, 


(B4) 
— 26.1053, 8.96339, —1.20016} . 
This formula is valid when 2 < R,/Ro < 193.7 and agrees 
with the Figure to within 7.996. 


C. DATA AND SOFTWARE AVAILABILITY 


The software and data required to reproduce our results are 
available under the digital object identifiers  10.5281/zen- 
0do.6368227 and 10.5281/zenodo.6371752, respectively. 
These repositories include the wind tunnel FLASH setup, 
the code we used to integrate the equation of motion of the 
engulfed companion, and the drag coefficients we measured 
in our hydrodynamical simulations. 
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